Method and Apparatus for Image Enhancement of Radiographic Images

ABSTRACT

A processing method for enhancing the image quality of an image, more particularly a digital medical grey scale image, that comprises the steps of a) decomposing an original image into multiple detail images at different resolution levels and/or orientations, b) processing the detail images to obtain processed detail images, c) computing a result image by applying a reconstruction algorithm to the processed detail ages, said reconstruction algorithm being such that if it were applied to the detail images without processing, then said original image or a close approximation thereof would be obtained, the processing of the detail images comprises the steps of: d) calculating at least one conjugate detail image, and e) computing at least one value of the processed detail images as a function of said conjugate detail image and said detail images.

TECHNICAL FIELD

The present invention relates to a method and a system for enhancing theimage quality of an image that is represented by a digital signal. Morein particular it relates to such method for use in a medicalradiographic imaging system, such as a digital radiography system.

Background of the Invention

In the ever continuing search for improved algorithms to enhance theimage quality of digital images, and especially of digital medicalimages, different techniques have been tried and tested. Theimprovements that are sought are obviously better enhancement results(such as better contrast enhancement, noise reduction, . . . ) orimproved computing performance with only limited artefacts ordistortions.

Commonly images represented by a digital signal such as medical imagesare subjected to image processing during or prior to displaying or hardcopy recording. The conversion of grey value pixels into values suitablefor reproduction or displaying may comprise a multi-scale imageprocessing method (also called multi-resolution image processing method)by means of which the contrast of the image is enhanced.

According to such a multi-scale image processing method an image,represented by an array of pixel values, is processed by applying thefollowing steps. First the original image is decomposed into a sequenceof detail images at multiple scales or resolution levels andoccasionally a residual image. Next, the pixel values of the detailimages are modified by applying to these pixel values at least oneamplification or conversion. Finally, a processed image is computed byapplying a reconstruction algorithm to the residual image and themodified detail images. Variations of this type of processing methodhave been described in the art, such as for instance in EP 1933273 B1.

There are limits for the behaviour of the amplification or conversionfunctions. Grey value transitions in the image can be distorted to anextent that the appearance becomes unnatural if the amplification orconversion functions are excessively non-linear. The distortions aremore pronounced in the vicinity of significant grey level transitions,which may result in overshoots at step edges and loss of homogeneity inregions of low variance facing strong step edges. The risk of creatingartefacts becomes more significant for CT images since they have sharpergrey level transitions, e.g. at the interface of soft tissue andcontrast media. Consequently, one has to be careful using themulti-scale techniques on CT images.

In European patent EP1933273 B1 a multi-scale contrast enhancementalgorithm which results in a contrast enhanced image while preservingthe shape of the edge transitions has been described. In one embodimentof the described method, translation difference images of at least oneapproximation image of the image are created at one or multiple scalesor resolution levels. Next, translation difference images arenon-linearly modified. Then at least one enhanced centre differenceimage at a specific scale is computed by combining modified translationdifference images at that scale or at a smaller scale.Spatially-localized phenomena derived from the image can be used tocreate enhanced center difference images. Finally, an enhanced image iscomputed by applying a reconstruction algorithm to the enhanced centerdifference images.

This method is, however, disadvantageous since non-linear modificationis applied to all translation difference images. Which might result inslow calculation.

In European patent EP2232434, the non-linear modification of the valuesof the translation difference image is steered by the values of anorientation map which comprises for each pixel on a specific scale alocal direction of interest or wherein values of weights that areapplied to the images that are summed up to form said amplificationimage are steered by the values of an orientation map which comprisesfor each pixel on a scale a local direction of interest.

This method is, however, disadvantageous since non-linear modificationis a function of only one local direction at a specific scale. Thismight be disadvantageous at corners or overlapping structures, oftenpresent in radiographic images.

Another multi-scale contrast enhancement algorithm which results in acontrast enhanced image while preserving the shape of the edgetransitions has been described in European application EP19209050.4. Inone embodiment of this method statistical measures of translationdifferences are used to enhance image. The last method is fast, but doesnot incorporate orientation in the conversion function to enhance theimage.

The application of the concept of steerable filters to image processingwas considered as a promising starting point to optimize imageprocessing algorithms with oriented filters. Oriented filters have beenfound useful in many vision and image processing tasks. One often needsto apply the same filter, rotated to different angles under adaptivecontrol, or wishes to calculate the filter response at variousorientations. A solution has been described in “The Design and Use ofSteerable Filters”, by Freeman and Adelson, IEEE transactions on patternanalysis and machine intelligence, Vol. 13, No. 9, September 1991, toprovide an efficient architecture to synthesize filters of arbitraryorientations from linear combinations of basis filters, allowing one toadaptively “steer” a filter to any orientation, and to determineanalytically the filter output as a function of orientation.

In addition, Freeman and Adelson use an approximation of Hilberttransform to form quadrature pairs. As such they can not only analysethe local orientation, but also the local amplitude (or magnitude) andphase.

Freeman and Adelson use the steerable filter to analyse images, andpropose following applications: analysing local orientation, angularlyadaptive filtering, contour detection and shape-from-shading analysis.They also introduce the concept of the steerable pyramid. Thecoefficients of the steerable pyramid image transform allow control overorientation analysis over all scales. Other pyramids based on quadraturemirror filters (QMF) or wavelets proved to be efficient for image codingapplications. The steerable pyramid is overcomplete, which limits itsefficiency but increase the convenience for image processing tasks. Theuse of steerable filters is optimal for orientation analysis. Thesimplest example of a steerable basis is a set of N (N−1)^(th)-orderdirectional derivatives of a circular symmetric function. In practice,according to theorem 3 of Adelson and Freeman, The Gaussian derivative,which can be written as the product of a Hermitian polynomial and aGaussian windowing function is suitable as a steerable filter. The setof basis filters are N rotations over n*π/N of such an (N−1)^(th) orderGaussian derivatives, with n∈{0, 1, . . . , N−1}. Any other rotation canbe computed as a linear combination of this set of basis filters.Adelson and Freeman also presented a set of separable steerable filtersas basis filters.

As explained in “The steerable pyramid: A flexible architecture formulti-scale derivative computation” by Simoncelli and Freeman, InProceedings International Conference on Image Processing (Vol. 3, pp.444-447). IEEE., October 1995, steerable filters are highly constrainedif they are used in a steerable pyramid. The steerable pyramid iscomputed recursively using convolution and decimation operations, and itis “self-inverting”. “Self-inverting” means that the matrixcorresponding to the inverse transformation is equal to the transpose ofthe forward transformation matrix, in the wavelet literature, such atransform is known as a tight frame. A major advantage of this steerablepyramid decomposition over other wavelet decompositions is that thesub-bands are translation- and rotation-invariant.

SUMMARY OF INVENTION

The invention provides a method for enhancing the contrast of anelectronic representation of an original image X represented by an arrayof pixel values by processing said image, as said out in independentclaim 1. The processing method of the invention comprises the steps ofa) decomposing said original image into multiple detail images atdifferent resolution levels and/or orientations, b) processing at leastone pixel of said detail images to obtain processed detail images, c)computing a result image by applying a reconstruction algorithm to theprocessed detail images, said reconstruction algorithm being such thatif it were applied to the detail images without processing, then saidoriginal image or a close approximation thereof would be obtained,characterized in that said processing of at least one pixel of saiddetail images comprises the steps of:

d) calculating at least one conjugate detail image, and e) computing atleast one value of the processed detail images as a function of saidconjugate detail image and said detail images.

The detail images can be calculated from the original image X, by firstcalculating a pyramid of increasingly lower-resolution images l₀, l₁, .. . , l_(K), each being a consecutively downscaled version of itspredecessor in said pyramid. Where the resolution of l₀ equals X andl_(K) is a low-pass residual image. And next by applying a convolutionto each of said lower-resolution images l_(k) ∈{l₀, l₁, l₂, . . . ,l_(K−1)} by a set of N steerable filters B_(n){B₀, B₁, . . . , B_(N−1)}to calculate a corresponding set of N detail images b_(kn) {b_(k0),b_(k1), . . . , b_(kN−1)} for each scale or resolution. B_(n) are Nrotated versions over n*π/N radians of filter B_(O). As put forward byFreeman and Adelson in the above cited article, a filter is suitable asa steerable filter on condition that any rotation of this filter can bewritten as a linear combination of the set of basis filters. Here, aconvolution of an arbitrary rotation of B₀ with l_(k) can be calculatedas a linear combination of b_(kn).

The steerable pyramid is computed recursively using convolution anddecimation operations, and it is “self-inverting”. A way to design suchfilters are described in “The steerable pyramid: A flexible architecturefor multi-scale derivative computation” by Simoncelli and Freeman, InProceedings International Conference on Image Processing (Vol. 3, pp.444-447). IEEE., October 1995. In addition, a high pass residual detailimage h_(O) is also calculated, because the detail images b_(kn) arebandpass filtered results. This can be calculated with filter H₀ or as aresidual of X-L₀.

For reconstruction, the inverse filters B′_(n) are applied to the detailimages. According to the self-inverting property, the inverse filtersare derived from the original filters by reflection of theircoefficients about the center. The advantage is that the steerablepyramid is self-inverting, and thus the errors introduced in thesubbands by quantization or non-linear modification, as proposed in thisinvention, will not appear as out-of-band distortions uponreconstruction. These filtered detail images are added with an upscaledversion of lower-resolution image I′_(k+1) resulting in an image I′_(k)which is, in case of unmodified detail images, a close approximation ofthe original l_(k).

The reconstructed image X′ is obtained by iteratively repeating theabove process for all scales k, starting from the low-pass residualimage l_(K), and finally adding the high pass residual h₀. X′ will be aclose approximation of X in case the detail images b_(kn) and h₀ areleft unmodified.

The steerable pyramid is an overcomplete representation. “Overcomplete”in this context means that the information that is present in the set ofresulting images is redundant. This reduces efficiency for instance forimage coding applications but it increases the convenience for imageprocessing tasks.

In addition, conjugate detail images c_(kn) are calculated byconvolution of said lower-resolution images l_(k) ∈{l₀, l₂, . . . ,l_(K−1)} with a set of N filters C_(n) {C₀, C₁, . . . , C_(N−1)}. WhereC₀ is the Hilbert transform of B₀ and where C_(n) is a rotation of C₀over n*π/N.

The invention differs from the state of the art in so far that a complexpyramid with detail images and conjugate detail images is constructedand that the detail images obtained from this complex steerable pyramidare modified before reconstruction. These modifications are non-linearand a function of b_(kn) and c_(kn) and can either selectively amplifyor suppress certain aspects of some pixels (i,j) in the detail images.After all detail images are processed, the reconstruction is performedbased on the inverse operation applied to all processed detail imagesb′_(kn) in order to obtain the processed and contrast enhanced versionof the grayscale image X′.

In practice one defines an amplitude A_(kn) and phase ϕ_(kn) for eachquadrature pair b_(kn) and c_(kn), and thus for each scale (k) andorientation (n):

$A_{kn} = \sqrt{b_{kn}^{2} + c_{kn}^{2}}$$\phi_{kn} = {a\tan( \frac{c_{kn}}{b_{kn}} )}$

The detail image can be defined as:

b _(kn)(i,j)=A _(kn)(i,j)*cos ϕ_(kn)(i,j)

According to the invention the amplitude is modified non-linearly into aprocessed amplitude A′_(kn), so that the processed detail images b′_(kn)are calculated by applying a modified amplitude and keeping the phase:

b′ _(kn)(i,j)=A′ _(kn)(i,j)*cos ϕ_(kn)(i,j)

With:

A′ _(kn)(i,j)=A _(kn)(i,j)g _(kn)(A ₀₀ . . . A _(0N−1) . . . A_(K−1N−1),ϕ₀₀ . . . ϕ_(0N−1) . . . ϕ_(K−1N−1))

and with g_(kn) being a positive amplification function.

It is the object of the present invention to provide a method forenhancing the image quality of an image, and more particularly indigital medical grayscale images. The present invention is advantageousin that it provides a method for enhancing the quality in a digitalimage so that the appearance remains natural in comparison to otheralgorithms where the result of the processing may become excessivelynon-linear, especially near sharp edges or lines, such that the methodreduces the distortions in the vicinity of significant grey leveltransitions.

Specific examples and preferred embodiments are set out in the dependentclaims.

The present invention is generally implemented as a computer programproduct adapted to carry out the method of any of the claims when run ona computer and is stored on a computer readable medium. The methods ofthe present invention can be applied for enhancing the image quality ofmedical images such as digital radiographic images, mammographic images,images obtained by computed tomography etc.

The embodiments of the systems and methods described herein may beimplemented in hardware or software, or a combination of both. However,preferably, these embodiments are implemented in computer programsexecuting on programmable computers each comprising at least one modulecomponent which comprises at least one processor (e.g. amicroprocessor), a data storage system (including volatile andnon-volatile memory and/or storage elements), at least one input device,and at least one output device. For example, and without limitation, theprogrammable computers (referred to here as computer systems) may be apersonal computer, laptop, personal data assistant, and cellulartelephone, smart-phone device, tablet computer, and/or wireless device.Program code is applied to input data to perform the functions describedherein and generate output information. The output information isapplied to one or more output devices, in known fashion.

Each program is preferably implemented in a high level procedural orobject oriented programming and/or scripting language to communicatewith a computer system. However, the programs can be implemented inassembly or machine language, if desired. In any case, the language maybe a compiled or interpreted language. Each such computer program ispreferably stored on a storage media or a device (e.g. ROM or magneticdiskette) readable by a general or special purpose programmablecomputer, for configuring and operating the computer when the storagemedia or device is read by the computer to perform the proceduresdescribed herein. The subject system may also be considered to beimplemented as a computer-readable storage medium, configured with acomputer program, where the storage medium so configured causes acomputer to operate in a specific and predefined manner to perform thefunctions described herein.

It is the object of the present invention to provide a method forenhancing the image quality of an image, and more particularly indigital medical grey scale images. Therefore, the coefficients of asteerable pyramid, i.e. the detail images, will be non-linearlymodified. The amplification or conversion function depends on scale,orientation, amplitude and phase.

The present invention is advantageous in that it provides a method forenhancing the quality in a digital image so that the appearance remainsnatural in comparison to other algorithms where the result of theprocessing may become excessively non-linear, especially near sharpedges or lines, such that the method according to the present inventionreduces the distortions in the vicinity of significant grey leveltransitions.

The present invention is further advantageous in that it provides amethod for reducing the noise level in a digital image so that excellentcontrast can be achieved without boosting the inherent image noise to adisturbing level.

Further advantages and embodiments of the present invention will becomeapparent from the following description and drawings.

BRIEF DESCRIPTION OF DRAWINGS

FIG. 1 depicts the steerable pyramid processing method by means of thedifferent filters and images involved. The original grayscale image X ishigh-pass filtered with filter H₀ into high-pass detail image h₀ and inlow-pass image l₀. The (in this case) three detail images b₀₀, b₀₁, b₀₂,can be calculated from the convolution of low-pass image l₀ with a setof (in this case) 3 steerable filters B_(n)∈{B₀, B₁, B₂}. Similarly,three conjugate images c₀₀, c₀₁, c₀₂, can be calculated from theconvolution of low-pass image l₀ with a set of (in this case) 3 filtersC, ∈{C₀, C₁, C₂}. Further low-pass filtering of l₀ and downscalingresults in l₁ from which in a similar way b₁₀, b₁₁, b₁₂, c₁₀, c₁₁, c₁₂and l₂ can be calculated. In the figure the process stops at residualimage l₂, but it can be further iterated and finally result in aresidual low-pass image l_(k) (l₂ in the example). The detail imagesb_(kn) are filtered with inverse filter B′_(n)∈{B′₀, B′₁, B′₂}. anditeratively added to the upscaled version of l_(k+1). In this way, andby adding high pass detail image h₀, the grayscale image X′ can bereconstructed.

FIG. 2 represents an example of filter sets that may be applied in themethod of the invention: low-pass filter L₀ and high-pass filter H₀ foruse at scale 0, low-pass filter L, 6 steerable filters B₀₋₅ rotated bysteps of 30 degrees and 6 conjugate filters C₀₋₅

FIG. 3 shows an example of detail images after applying filters fromFIG. 2 .

FIG. 4 shows an example of amplitude images by using filters from FIG. 2.

DESCRIPTION OF EMBODIMENTS

In the following detailed description, reference is made in sufficientdetail to the above referenced drawings, allowing those skilled in theart to practice the embodiments explained below.

In the following, a more detailed overview of the preferred embodimentswill be explained in more detail. As a first step in the process(reference is made to FIG. 1 with a complete overview of the process),and before the decomposition steps are started, a high-pass filter H₀ isapplied to the original image X, and results in the high-pass filtereddetail image h₀. Since detail images b_(kn) are band-pass filtered, thisdetail image h₀ is needed to reconstruct the finest detail. In apreferred embodiment, H₀ is chosen in such a way that it isself-inverting, and has to be applied again before reconstruction.

The high-pass filtering step is performed as follows: a) a high-passfilter operator is applied to the original image X, obtaining ahigh-pass filtered detail image h₀ and storing said image h₀ in memory,b) a low-pass filter operator L₀ is applied to the original image X,obtaining a low-pass filtered image l₀ at the original image resolutionand storing said low-pass filtered full-resolution image l₀ in a memory.

In another embodiment, the high-pass detail image h₀ could be obtainedas the difference between X and low-pass image l₀.

As a second step, the low-pass filtered full-resolution image l₀ of theoriginal image X is scaled down into a pyramid of increasinglylower-resolution images l_(I), l₂, . . . , l_(K), each consecutive imagebeing a low-pass filtered (with filter L) and downscaled version of itspredecessor in said pyramid. By downscaling the original size image l₀,one thus obtains a smaller sized next image l₁ that may be furtherscaled down.

The downscaling of the low-pass filtered image l_(k) can be performed byany known downscaling algorithm such as bi-linear interpolation,box-scaling, bi-cubic interpolation or alike. A preferred downscalingmethod for the application in the invention would be one where thespatial resolution is halved along both axes of the image. The amount ofpixels of the remaining image would therefore be reduced by 4.

Since a downscaling operation on a digital image can be considered as asmoothing operation, it means that the processing of a downscaled imageis equivalent to the processing of a smoothed version of the image. Inthis respect, it can be understood that the pyramidal decompositionapproach leads to the image processing of the same image at differentresolution levels or applying to different signal frequencies present inthe image.

In a preferred embodiment the filters depicted in FIG. 2 and proposed bySimoncelli and Freeman, combined with reduction of resolution by 2×2 areused.

In a third step, and after the decomposition of the original image intolower resolution images l₀, l₁, l₂ . . . l_(k), each resolution versionl_(k)∈{l₀, l₁, l₂ . . . l_(K−1)} is further decomposed by applying a setof N filters B_(n) ∈{B₀, B₁, . . . , B_(N−1)}, where B_(n) is a rotatedversion of B₀ by a rotation of n/N*π.

This is achieved by applying a convolution of this set of filters toeach of the lower resolution images l_(k). As such for each l_(k), Ndetail images b_(kn) are obtained and each of the N filters will gaugethe local orientation at lower resolution image l_(k). As a result, foreach detail image b_(kn), k indicates the scale or resolution level, andn the orientation.

In the preferred embodiment, this set of filters are steerable. Thisimplies that, the convolution of l_(k) with any rotation of B₀ can becomputed as a linear combination of basis images b_(kn).

The simplest example of a steerable basis is a set of N (N−1)^(th)-orderdirectional derivatives of a circularly symmetric function, for examplethe Gaussian derivatives that are represented by a square pixel kernelof a manageable size. Since the kernels need to be convolved with theimage, it speaks for itself that the choice of the kernel size will be atrade-off between accuracy and computing speed. A larger convolutionkernel will require more calculation time, but will render more accurateresults. Typical higher order Gaussian derivatives need larger kernelsizes. A typical kernel size will be between 5×5, 7×7, 9×9 or 11×11.

The steerable pyramid is computed iteratively using convolution anddecimation operations, and it is “self-inverting”. A way to design suchfilters is described in “The steerable pyramid: A flexible architecturefor multi-scale derivative computation” by Simoncelli and Freeman, InProceedings International Conference on Image Processing (Vol. 3, pp.444-447). IEEE., October 1995. These filters need to be band-limited,have a flat system response and can be used iteratively in order toconstruct a pyramid. In addition, they need to be steerable andself-inverting.

The derivative order (N−1) can be chosen. Higher order derivatives willresult in better orientation resolution but also in larger kernels andmore (N) detail images and will consequently require a highercomputational effort. For instance, a first order directional derivativewill not be able to discriminate between more than two orientations(only steps of 90 degrees), and is thus not useful for enhancingstructures with multiple orientations. Since multiple orientations arecommonly present in radiographic image, because of its additive nature,higher order derivates are preferred.

In practice, and in a preferred embodiment, a steerable filter-sethaving 6 orientations with steps of 30 degree rotation are taken, tomatch the simple cells in the primary visual cortex of the human visualsystem. An illustration of such a steerable filter set is given in FIG.2 as B₀₋₅.

The convolution operation in the example of each of these filters withthe lower resolution image l_(k) results in detail images b_(k0-5) (anequal amount as the number of steerable filters in the set), of which anexample is shown in FIG. 3 . For each convolution with filters B₀₋₅ acorresponding detail image b_(k0-5) is produced. In one embodiment thefilters B_(n) are calculated by a rotation of filter B₀ with n/N*π.

In another embodiment separable steerable filters are used to calculatethe rotated version or a close approximation. Convolution with separablebasis filters is computational cheaper. For most steerable filters, thebasis filters are not x-y separable. According to Freeman and Adelson,Appendix D, for each filter that can be written as polynomial in x andy, there is a x-y separable basis. These basis functions can then beused to steer the filter B₀ to any wanted rotation, and more specific torotations n/N*π.

At the same time, or in a fourth step, the Hilbert transform of theabove steerable filter-set B_(n) (indicated as C₀₋₅ in the example ofFIG. 2 ) is convolved with each lower resolution image l_(k) tocalculate the conjugate detail images C_(k0-5). These conjugate detailimages comprise the “imaginary” part of the signal in the image.

In one embodiment, C₀ is calculated as Hilbert transform of B₀, andC_(n) are rotated versions over n/N*π of C₀. In another embodiment, C₀is a polynomial approximation of the Hilbert transform of B₀. A set ofN+1 steerable filters is created as described by Adelson and Freeman. Inthis approach, the filter response that corresponds to rotation of n/N*πis calculated as a linear combination of a set of (separable) basisfilters.

The Hilbert transform is important in signal processing to find theanalytic representation of a real-valued signal u(t), or the so-calledanalytic signal. Specifically, the Hilbert transform of u(t) is itsharmonic conjugate v(t), a function of the real variable tsuch that thecomplex-valued function u(t)+iv(t) admits an extension to the complexupper half-plane satisfying the Cauchy—Riemann equations.

In practice, by calculating the Hilbert transform of the real part of adigital signal, one can find the (complex valued) analytical signal. Theamplitude of this analytical signal corresponds to the ‘energy’ of thesignal because it is not modulated. More specifically, by applying aRiesz transform, which is an N-dimensional generalization of the Hilberttransform, it is possible to find the amplitude, phase and orientationof the signal at a certain scale (or “level”) of the image (and not onlythe real response). Calculating these Riesz transforms at differentscales of the image will provide knowledge about the structure in theimage at the different scales, which can then be used to modify thedetail images accordingly. In this invention the detail images, or thereal part of the steerable pyramid, will be modified, based on thecalculated amplitude, phase and orientation information, and from whichan enhanced image X′ can be calculated.

In order to explain this further, a simple example is given to show theimportance of phase on image enhancement. For an odd order set offilters B_(n), the “real” part of the signal, or the detail imageb_(kn), will generate more signal for edges than for lines; the detailimages b_(kn) thus will mainly comprise signals representing thepresence of edges in the lower resolution image l_(k). For odd orderB_(n), the Hilbert transform C_(n) has even order, and the “imaginary”part of the signal, at phase +90°, c_(kn) will generate more signal forlines than for edges, and the conjugate detail image c_(kn) will mainlycomprise signals representing lines. In this way, the phase of oursignal will be a measure for the “edgeness” (i.e. the amount of edgespresent in a pixel) in the lower resolution image. And the amplitude ormagnitude of the signal gives a measure for the “energy” in the signal,regardless of whether this signal represents a line or edge. Inaddition, one can obtain this amplitude and phase for each orientationand scale. As such, one can control the modification of detail imageb_(kn) based on the amplitude and ‘line-edge-ness’ (phase) in eachdirection, and at finer scales one could for instance enhance edges morethan lines because the former are less prone to noise. In a simple andpreferred embodiment, the amount of ‘line-edge-ness’ or phaseenhancement is rotationally invariant or independent of its direction.

In a preferred embodiment, amplitude A_(kn) and phase ϕ_(jai) for eachquadrature pair b_(kn) and c_(kn), and thus for each scale (k) andorientation (n) are calculated as:

$A_{kn} = \sqrt{b_{kn}^{2} + c_{kn}^{2}}$$\phi_{kn} = {a\tan( \frac{c_{kn}}{b_{kn}} )}$

Where the amplitude A_(kn) is a measure for contrast at a certain scalek and orientation n, while the phase ϕ_(kn) is a measure forline/edgeness at a certain scale k and orientation n.

The corresponding detail image is equal to:

b _(kn)(i,j)=A _(kn)(i,j)*cos ϕ_(kn)(i,j)

It is the object of this invention to non-linearly modify the amplitudeinto a processed amplitude A′_(kn), in order to obtain processed detailimages with the same phase but a modified amplitude (contrast):

b′ _(kn)(i,j)=A′ _(kn)(i,j)*cos ϕ_(kn)(i,j)

wherein:

A′ _(kn)(i,j)=ƒ_(kn)(A ₀₀ . . . A _(0N−1) . . . A _(K−1N−1),ϕ₀₀ . . .ϕ_(0N−1) . . . ϕ_(K−1N−1))

f_(kn) being a conversion function. Typically, ƒ_(kn) is a non-linearfunction in order to obtain contrast enhancement in the reconstructedimage.

Further is it more convenient to define this conversion function as amultiplicative amplification function or a gain function g_(kn):

A′ _(kn)(i,j)=A _(kn) g _(kn)(A ₀₀ . . . A _(0N−1) . . . A _(K−1N−1),ϕ₀₀. . . ϕ_(0N−1) . . . ϕ_(K−1N−1))

In a preferred embodiment this amplification function is translation-and rotation-invariant (such as is the steerable pyramid itself). Inother words, translation or rotation of image X has almost no impact andthe processing result. Or in other words, any rotation or translation ofimage X followed by decomposition, modification of the detail images,reconstruction and inverse rotation or translation should have no oralmost no effect on enhanced image X′. In another embodiment, theamplification function is translationally invariant and onlyquasi-rotationally invariant. In this case, there is no or almost noeffect on the processing for some specific rotations, for instance stepsof 30 degrees.

In a preferred embodiment, this amplification is only depending on thephase and amplitude at the same scale k:

g_(kn)(A₀₀…. A_(0N − 1)… A_(K − 1N − 1), ϕ₀₀…. ϕ_(0N − 1) …ϕ_(K − 1N − 1)) = g_(kn)(A_(k0)….A_(kN − 1), ϕ_(k0)….ϕ_(kN − 1))

In another embodiment, g_(kn) is not only depending on amplitude andphase of the scale k. It could also be a function of lower scales. E.g.one could calculate correlation between A_(kn) and A_(k+1n). Pixels atdifferent scale with high correlations could get similar amplificationto reduce artefacts.

In practice, this amplification function is split in two or moremultiplicative amplification functions. In a preferred embodiment g_(kn)is split in three parts:

g _(kn)(A _(k0) . . . A _(kN−1),ϕ_(k0) . . . ϕ_(kN−1))=g1_(kn)(A_(kn))g2_(kn)(ϕ_(kn))g3_(kn)(A _(k0) . . . A _(kN−1))

Where g1_(kn) is amplitude, g2_(kn) is phase and g3_(kn) is orientationdependent.

In one specific embodiment the first amplitude dependent amplificationfunction can be defined as:

${g1_{kn}} = \frac{f1_{kn}( A_{kn} )}{A_{kn}}$

Wherein A_(kn) is the amplitude for scale k and orientation n and is ameasure for contrast, and mapping function ƒ1_(kn) can be defined as forinstance illustrated in EP 527 525, p 9 line 35:

ƒ1_(kn) :x→αx ^(p)

In this embodiment x will be A_(kn). (rather than being defined as thedetail image pixel as in EP 527 525). As such:

g _(kn) =αA _(kn) _(p) ^(p−1)

where the power p is chosen within the interval 0<p<1, preferably withinthe interval 0.5<p<0.9. A comparative evaluation of a large number ofcomputed radiography images of mainly thorax and bony structures by ateam of radiologists indicated that p=0.7 is the optimal value in mostcases. α specifies a gain factor.

If this multiplicative amplification factor is applied to a detail image

b _(kn)(i,j)=g _(kn)(A _(kn)(i,j))*A _(kn)*cos ϕ_(kn)(i,j)

then all details with a low amplitude will be boosted relative to theimage details which originally have a higher amplitude.

In this respect, the above power function ƒ1_(kn): x→αx^(p) proved toperform very well, but it is clear that an infinite variety ofmonotonically increasing mapping functions can be found that willenhance subtle details with low amplitude. The main requirement is thatg_(kn) is higher, and thus the slope of said mapping function issteeper, in the region of argument values that correspond to lowelementary contrast.

In an alternative embodiment, excessive noise amplification can beavoided by using a composite mapping function:

$ {f_{kn}:x}arrow{{\alpha{c^{p_{2}}( \frac{x}{c} )}^{p_{1}}{if}x} < c} $x → αx^(p₂)ifx ≥ c

where the power p₂ is chosen in the interval 0<p₂<1, preferably0.5<p₂<0.9, and most preferably p₂=0.7 (however the preferred value ofp₂ highly depends upon the exact kind of radiological examination, andmay vary). The power p₁ in this equation should not be smaller than p₂:p₁≥p₂, where the cross-over abscissa c specifies the transition pointbetween both power functions.

Decreasing the power p₂ will further enhance the contrast of subtledetails, but at the same time the noise component will also beamplified. The noise amplification can be limited by choosing a powervalue p₁ larger than p₂, preferably 1.0, so that the slope of themapping function is not extremely steep for the range of very smallabscissae.

Ideally the cross-over abscissa c should be proportional to the standarddeviation of the noise component (assuming additive noise), so thelowest amplitude details buried within the noise along with the majorityof the noise signals will only be moderately amplified, according to thevalue of the functional part controlled by the power p₁, while thedetail signals just exceeding the noise level will be amplified muchmore according to the value of the functional part controlled by thepower p₂. The decreasing value of the latter functional part stillensures that the subtle details above the noise level are boostedrelative to the high amplitude details. In this respect, the abovecomposite power function proved to perform very well, but it is clearthat an infinite variety of monotonically increasing mapping functionscan be found that will enhance subtle details without boosting the noiseto an excessive level.

The main requirement is that g_(kn) is higher, and thus the slope ofsaid mapping function is steeper, in the region of argument values thatcorrespond to low elementary contrast than it is either in the sub-rangeof very small elementary contrast which mostly correspond to noise, orin the range of the larger elementary contrast.

When all detail images of the decomposition are modified using the samemapping according to one of the above methods, a uniform enhancementover all orientations, and a rotation invariant conversion is obtained.In addition, the enhancement will also be uniform over all scales. In aslightly modified embodiment, where a different mapping function is usedat each resolution level e.g. by multiplying one of the above describedmapping functions or multiplicative factor with a resolutionlevel-dependent coefficient, it is possible to further increase thesharpness by setting the coefficient corresponding to the finestresolution level to a substantially higher value than the othercoefficients as follows:

y=A _(i) F(x)

for i=0 . . . K−1

where F(x) is one of the above described mappings, K is the number ofresolution levels, and A_(i) is a level dependent coefficient, e.g.A₀>1, and A_(i)<1 for 1≤i≤K−1.

In one specific embodiment the second, phase dependent, amplificationfunction can be defined as:

${g2_{kn}( {\phi_{kn}( {i,j} )} )} = {\frac{1}{n_{s}}( {1 - {{\phi_{kn}( {i,j} )}*\frac{2}{\pi}}} )^{2}}$

In which n_(s) is a noise suppression factor. A higher value of n_(s)results in more attenuation, a lower value in more amplification.g2_(kn) may also be limited to a band of values defined between an upperand lower bound. Especially for those functions with high values.Typical values for lower and upper bound are found to give good resultswhen they are respectively 0.5 and 1. A smaller value of the lower bound(e.g. zero) will result in more attenuation of noise and cartooning likeartefacts, while larger values of the upper bound (e.g. 2) will resultin increased enhancement of the edges, and overshoot artefacts at theedges. n_(s) is chosen for each scale. For finer scales, with typicallymore noise, a larger noise suppression parameter is chosen, for coarserscales a smaller noise suppression is chosen. As an example, the noisesuppression parameter n_(s)=1.2 at scale 0, whereas n_(s)=0.5 at scale 1results in natural looking noise suppression in the reconstructed image.

The role of this function is to attenuate detail image pixelsrepresenting lines or points, because they are probably, especially atfiner scales, noise. Lines have a phase close to 90 degrees, edges at 0and 180 degrees for odd-ordered B_(n) filters. The function g2_(kn)( )needs to have a lower amplification at a phase of 90 degrees than at 0degrees. Since it is a multiplicative amplification function, it needsto be positive. In addition, this function also needs to be symmetricaround 90 degrees, e.g same amplification has to be found for 90−x asfor 90+x. Note that for even-ordered B_(n) filters these functions needto be shifted with 90 degrees. The above function proved to perform verywell, but it is clear that an infinite variety of functions can be foundthat meet these requirements.

Some examples of such functions are given below:

${g2_{kn}( {\phi_{kn}( {i,j} )} )} = {\frac{1}{n_{s}}\frac{1}{\sin( {\phi_{kn}( {i,j} )} )}}$${g2_{kn}( {\phi_{kn}( {i,j} )} )} = {\frac{1}{n_{s}}\frac{1}{{\tan( {\phi_{kn}( {i,j} )} )}^{2}}}$${g2_{kn}( {\phi_{kn}( {i,j} )} )} = {\frac{1}{n_{s}}{\cos( {\phi_{kn}( {i,j} )} )}^{2}}$

In another embodiment, n_(s), as well as the lower and upper bound maybe calculated based on particular parameters that can be derived fromthe image data (or detail image data). For instance, histogram values ofb_(kn)(i,j) and c_(kn)(i,j) or a combination of such parameters may beapplied in a mathematical function to derive the parameters above(n_(s), and lower and upper bound).

In one specific embodiment the third, orientation dependent,amplification function can be defined as a function:

$\begin{matrix}{{g3_{kn}( {A_{k0}\ldots A_{{kN} - 1}} )} = {N( {{{if}n}=={\arg\max\limits_{n}A_{kn}}} )}} \\{= {0({elsewhere})}}\end{matrix}$

In this way, A′_(kn) will be only non-zero for orientation n with thelargest amplitude. This will reduce noise but tends to introduce manyartifacts into the image. Especially at corners and overlappingstructures. This could be improved in a specific embodiment by using:

$\begin{matrix}{{g3_{kn}( {A_{k0}\ldots A_{{kN} - 1}} )} = {( {1 - \alpha} ) + {N{\alpha( {{{if}n}=={\arg\max\limits_{n}A_{kn}}} )}}}} \\{= {( {1 - \alpha} )({elsewhere})}}\end{matrix}$

In this function the direction with maximum amplitude

$\underset{n}{\arg\max}A_{kn}$

is calculated for each pixel, and g3_(kn) is larger for n close

$\underset{n}{\arg\max}{A_{kn}.}$

In addition, sum of amplification over all directions is kept the same.

With α is a parameter. If α=0 than g3_(kii) is identity function, if α=1only orientation with maximum amplitude are kept (i.e. identical toprevious embodiment). α controls the amount of ‘orientation smoothing’,higher α will give more smoothing but more artifacts. These artifactsare more prominent at corners or overlapping structures. α can bedifferent for different scales, typical for finer scales which containmore noise, a higher α is chosen.

In addition, this function is not completely rotationally invariant, butquasi rotationally invariant (only for rotations of n/N*π). In animproved embodiment, this could be solved by interpolation. In aspecific embodiment one could find n_(max) by argmax_(n) A_(kn) withsub-integer accuracy by using e.g. quadratic interpolation. In anotherembodiment, n_(max) is found with sub-integer accuracy by calculatingamplitudes for more rotations than only n*π/N.

$\begin{matrix}{{g3_{kn}( {A_{k0}\ldots A_{{kN} - 1}} )} = {( {1 - \alpha} ) + {N\alpha*( {{{ceil}( n_{\max} )} - n} )( {{{if}n}=={{floor}( n_{\max} )}} )}}} \\{= {( {1 - \alpha} ) + {N\alpha*( {n - {{floor}( n_{\max} )}} )( {{{if}n}=={{ceil}( n_{\max} )}} )}}} \\{= {( {1 - \alpha} )({elsewhere})}}\end{matrix}$

In another embodiment, the steering equation, as for instance given inFreeman and Adelson, can be used instead of interpolation. E.g. for5^(th) order polynomial basis filter:

g3_(kn)(A_(k0)….A_(kN − 1)) = (1 − α) + N * α * k_(n)(θ)${{With}:{k_{n}(\theta)}} = {\frac{1}{6}( {{2*\cos( {\theta - \theta_{n}} )} + {2*\cos( {3*( {\theta - \theta_{n}} )} )} + {2*\cos( {5*( {\theta - \theta_{n}} )} )}} )}$${{With}:\theta} = {n_{\max}*\frac{\pi}{N}}$${{and}{}\theta_{n}} = {n*\frac{\pi}{N}}$

It is clear that many functions can be constructed wherein amplitudes inorientations n close to the maximum amplitude orientation n_(max) aremore enhanced than other orientations.

In another embodiment, multiple local maximum directions could be foundand enhanced. This could be advantageous at locations with multipleorientations, such as crossing structures.

In another embodiment, one could make α a function of the spread of mainorientations of A_(kn) (=argmax_(n)A_(kn)) in a M×M neighborhood aroundpixel (i,j). The goal is to reduce α at corners and overlappingstructures, and thus reduce artefacts. These points typically haveneighboring points with many orientations and thus a large spread.Points at clearly oriented lines and edges have a much smaller spread.In one specific and preferred embodiment one could calculate the mean ofall orientations in a M×M neighborhood:

$R = \sqrt{( {\sum\limits_{i \in {M \times M}}( {w_{i}\sin\frac{2\pi}{N}\arg\max_{n}A_{kn}} )} )^{2} + ( {\sum\limits_{i \in {M \times M}}( {w_{i}\cos\frac{2\pi}{N}\arg\max_{n}A_{kn}} )} )^{2}}$

With w_(i) is a weight so Σ_(i∈M×M)(w_(i))=1 (e.g. for 3×3: w_(i)=1/9)

R equals 1 (R=1) in the case that angles

$\frac{2\pi}{N}\arg\max_{n}A_{kn}$

are the same Tor all pixels in M×M. R equals 0 (R=0) in the case thatangles are distributed evenly along unit circle.

One could define circular standard deviation stdev as:

stdev=√{square root over (−2*log(R))}

In a specific and preferred embodiment, one could make a dependent onthis stdev. If stdev is close to zero, one wants to get maximal focusand thus a large α, otherwise a smaller α is preferred. It is clear thatmany variations exist with similar functionality, in a specific andpreferred embodiment:

$\alpha = {\alpha_{0}\exp( {- \frac{stdev^{2}}{2*\sigma_{0}^{2}}} )}$

wherein suitable values are:

$\sigma_{0} = \frac{\pi}{N}$

α₀=0.5 at scale 0, 0.15 at other scales.

In another embodiment, α or α₀ depends on the gray value of the originalimage or on a segmentation mask. In medical radiographic images we foundthat it is advantageous to take different α₀ for different tissue types,e.g. a larger α₀ for bone. This could be defined by the original greyvalues, but this is often not a good approach since greyscale values maybe similar for different tissue types.

Defining this α₀ based on the final or an approximation of the finalpixel value is often a better approach. This could be done during asecond processing or by reconstruction of scale k+1. Another approachcould be to use a segmentation map. This segmentation could be forinstance a segmentation map of different tissue types, obtained by forinstance a convolutional neural network.

In another specific embodiment one could train or fit a regression for:

g _(kn)(A ₀₀ . . . A _(0N−1) . . . A _(K−1N−1),ϕ₀₀ . . . ϕ_(0N−1) . . .ϕ_(K−1N−1))

One could apply this multiplicative amplification function incombination with other multiplicative amplification functions. In onespecific embodiment, g_(kn) is a combination of amplitude dependentamplification function g1_(kn)(A_(kn)) and a trained amplificationfunction gREG_(kn). In this way one could use g1_(kn) to modify theamplitude and thus the contrasts (as explained before), and train afunction for gREG_(kn) to reduce artifacts and/or noise:

g_(kn)(A₀₀…. A_(0N − 1)… A_(K − 1N − 1), ϕ₀₀…. ϕ_(0N − 1)… ϕ_(K − 1N − 1)) = g1_(kn)(A_(kn)) * gREG_(kn)(A₀₀….A_(0N − 1)…A_(K − 1N − 1), ϕ₀₀….ϕ_(0N − 1)…ϕ_(K − 1N − 1))

In another embodiment, gREG_(kn) is applied to reduce noise and/orartifacts in a processed amplitude (gREG_(kn)*A_(kn)), next thisprocessed amplitude is used in another amplification function, e.g.g1_(kn).

g_(kn)(A₀₀…. A_(0N − 1)… A_(K − 1N − 1), ϕ₀₀…. ϕ_(0N − 1)… ϕ_(K − 1N − 1)) = g1_(kn)(gREG_(kn) * A_(kn)) * gREG_(kn)(A₀₀….A_(0N − 1)…A_(K − 1N − 1), ϕ₀₀….ϕ_(0N − 1)…ϕ_(K − 1N − 1))

In order to train this regression gREG_(kn), a training set of originalimages X(i,j) is collected to which artificial noise or artefacts areadded to obtain a set of artificial images X_(n)(i,j). This added noisemay be Poisson noise, Gaussian noise, other noise or any combination.Artefacts could be dead pixels/lines, a simulated overexposure, gridlines, ringing artefacts, scatter or other artefacts or any combinationthereof. The original images X(i,j) may be a collection of high doseimages, phantom images or a lower resolution images with less noise. Foreach X(i,j) A_(kn) and for each X_(n)(i,j) An_(kn) is calculated, next

${gREG_{kn}} = \frac{A_{kn}}{An_{kn}}$

is calculated. Next a set of amplitudes and phases (An₀₀ . . . An_(0N−1). . . An_(K−1N−1), ϕn₀₀ . . . ϕn_(0N−1) . . . ϕn_(K−1N−1)) are computedfor each X_(n)(i,j) as explained before, and the wanted amplificationgREG_(kn) is fitted for each pixel to this set of amplitudes and phases.The goal is to train or fit a function fin such a way that:

f:(A ₀₀ . . . A _(0N−1) . . . A _(K−1N−1),ϕ₀₀ . . . ϕ_(0N−1) . . .ϕ_(K−1N−1))→gREG_(kn)(gREG_(kn) *An _(kn) −A _(kn))² is minimized

In an another embodiment, one could not only take the amplitudes andphases for pixel (i,j), but also from other pixels. In a specificembodiment, the pixels in a M×M neighbourhood are taken into account.

In order to train f one can use linear regression, Bayesian networksprincipal component, support vector machine, neural network, or othermachine learning algorithms. This function can then be applied on realdata to reduce noise/artefacts.

High-pass filtered detail image h₀ may be transformed into an enhancedhigh-pass filtered detail image h′₀, by means of existing methods basedon the pixel value, surrounding pixel values or by calculatingstatistical measures as further described in for instance EP1933273 B1,in European applications EP19202349.7 and EP19209050.4. In anotherembodiment the values of (b_(kn), c_(kn)) can be used to modify thishigh pass filtered detail image h₀. For example, for each pixel, themean of amplifications at scale 0 (g_(0n))

$( {\sum_{n \in N}\frac{g_{on}}{N}} )$

can be taken to amplify detail image pixel values h₀(i,j). In anotherembodiment, a hybrid method is used, based on existing methods andvalues of (b_(kn), c_(kn)), e.g. amplification can be defined as aweighted average of both techniques. In another embodiment, theprocessed detail image h′₀ will not be calculated and will be left outfor reconstruction, because h₀ may mainly contain noise.

Low pass residual l_(K) may be transformed into an enhanced low-passresidual I′_(K) by means of state of the art. One could for instance useother techniques as for instance described in EP1933273 B1, in Europeanapplications EP19202349.7 and EP19209050.4. In a preferred embodiment,the detail images b_(kn) are only calculated for the first scales,because orientation is not very important for lower scales. The low-passresidual I′_(K) is then enhanced by for instance a multi-scale contrastenhancement algorithm based on statistical measures of translationdifferences as described in European application EP19209050.4. Inanother embodiment, the low-pass residual will not be processed at all(identity processing) or will be left out for the reconstruction.

After processing the detail images b_(kn), the processed detail imagesare obtained as b′_(kn). Subsequently, the processed detail images areconvolved with the inverse filter set B′_(n). Because the steerablepyramid is self-inverting, these inverse filters are derived from theoriginal filters B_(n) by reflection of their coefficients about thecenter. In order to reconstruct the processed image X′, all processeddetail images are summed up (with filtering), including the residuallow-pass filtered version I′_(K), and the high-pass filtered version h′₀(with filtering).

$X^{\prime} = {l_{K}^{\prime} + {\sum\limits_{k = 0}^{K - 1}{\sum\limits_{n = 0}^{N - 1}{B_{n}^{\prime} \otimes b_{kn}^{\prime}}}} + {H_{0} \otimes h_{0}^{\prime}}}$

The detail images are modified starting from the lowest resolutiondetail image up to the finest level. Subsequently, the reconstructionalgorithm is applied to the modified detail images b′_(kn), againstarting with the lowest resolution detail image up to the finest level.

In addition, other enhancements may be applied to b_(kn) prior to orafter the application of the conversion operator g_(kn)(i,j).

This could be generalized to more dimensions and especially to 3D imagesX as proposed by Freeman and Adelson. This can be advantageous in theprocessing of temporal data (e.g. fluoroscopic radiographic imaging) or3D data sets such as CT, cone beam CT, digital tomosynthesis and MRI. 3Dfilters which can be written as N^(th) order polynomials multiplied witha radial windowing function (e.g. Gaussian) can be steered by a set of(N+1)² basis filters. In order to keep computational effort low, the useof x-y-z separable filters is even more advantageous. Similar toSimoncelli and Freeman such set of steerable filters that are used in asteerable pyramid can be computed. The Hilbert transform is used tocalculate conjugate detail images. Detail image and conjugate detailimages are used to calculate amplitude and phase for each orientationand nonlinear modification is applied before reconstruction.

1-24. (canceled)
 25. A method for enhancing the contrast of a digitalrepresentation of an original image X represented by an array of pixelvalues by processing said image, said processing comprising the stepsof: a) decomposing said original image X into a high-pass filtereddetail image h₀, a set of detail images b_(kn) at multiple resolutionlevels (K) and/or multiple orientations (N), b) processing at least onepixel (i,j) of said detail images h₀ and b_(kn) to obtain processeddetail images h′₀ and b′_(kn) c) computing a result image X′ by applyinga reconstruction algorithm to said processed detail images h′₀ andb′_(kn), d) calculating at least one conjugate detail image c_(kn), e)computing at least one pixel value of the processed detail images h′₀and b′_(kn) as a function of said conjugate detail image c_(kn) and saiddetail images h₀ and/or b_(kn) wherein a pixel of said processed detailimages is calculated as:b′ _(kn)(i,j)=A′ _(kn)(i,j)*cos ϕ_(kn)(i,j) wherein ϕ_(kn) is a phasedetail image, calculated as$\phi_{kn} = {a\tan( \frac{c_{kn}}{b_{kn}} )}$ and, whereinA′_(kn) is an amplitude of said processed detail image, calculated asA′ _(kn)(i,j)=ƒ_(kn)(A ₀₀ . . . A _(0N−1) . . . A _(K−1N−1),ϕ₀₀ . . .ϕ_(0N−1) . . . ϕ_(K−1N−1)) wherein, {A₀₀ . . . A_(0N−1) . . .A_(K−1N−1)} is an amplitude of at least one pixel of detail images, and{ϕ₀₀ . . . ϕ_(0N−1) . . . ϕ_(K−1N−1)} is a phase of at least one pixelof detail images, and wherein the amplitude of a detail image A_(kn){A_(k0), A_(k1), . . . , A_(kN−1)} is calculated as:A _(kn)=√{square root over (b _(kn) ² +c _(kn) ²)}.
 26. The method ofclaim 25, wherein detail images b_(kn) are calculated by: a) downscalingthe original image X into a pyramid of increasingly lower-resolutionimages l₀, l₁, . . . , l_(K−1), each being a consecutively downscaledversion of its predecessor in said pyramid, and b) applying aconvolution to each of said images l_(k) ∈{l₀, l₁, l₂, . . . , l_(K−1)}by a set of N filters B_(n) {B₀, B₁, . . . , B_(N−1)} to calculate acorresponding set of N detail images b_(kn) {b_(k0), b_(k1), . . . ,b_(kN−1)}, wherein each B_(n) is a rotated version of filter B_(o). 27.The method of claim 26, wherein said N filters B_(n) are steerablefilters, in such a way that a convolution of an arbitrary rotation of B₀with said lower-resolution image l_(k) can be calculated as a linearcombination of said detail images b_(kn).
 28. The method of claim 26,where said N filters B_(n) are calculated as a linear combination ofsteerable filters.
 29. The method of claim 26, wherein N is differentfor each resolution level.
 30. The method of claim 25, wherein aconjugate detail image c_(kn) is calculated by downscaling the originalimage X into a pyramid of increasingly lower-resolution images l₀, l₁, .. . , l_(K−1), each being a consecutively downscaled version of itspredecessor in said pyramid and applying a convolution to each of saidimages l_(K)∈{l₀, l₁, l₂, . . . , l_(K−1)} by a set of N conjugatefilters C₀, C₁, . . . , C_(N−1).
 31. The method of claim 30, wherein Nconjugate filters C₀, C₁, . . . , C_(N−1) are rotations of a Hilberttransform of said filter B₀.
 32. The method of claim 30, wherein Nconjugate filters C₀, C₁, . . . , C_(N−1) are calculated as a linearcombination of steerable filters.
 33. The method of claim 25, whereinf_(kn) is rotationally invariant.
 34. The method of claim 25, whereinA_(kn)^(′)(i, j) = f_(kn)(A₀₀….A_(0N − 1)…A_(K − 1N − 1), ϕ₀₀…. ϕ_(0N − 1) …ϕ_(K − 1N − 1)) = A_(kn) * g_(kn)(A₀₀….A_(0N − 1)…A_(K − 1N − 1), ϕ₀₀…. ϕ_(0N − 1) …ϕ_(K − 1N − 1))and wherein g_(kn) is multiplicative positive amplification function.35. The method of claim 34, wherein:A′ _(kn)(i,j)=A _(kn)(i,j)g1_(kn)(A _(kn))g2_(kn)(ϕ_(kn))g3_(kn)(A ₀₀ .. . A _(0N−1) . . . A _(K−1N−1),ϕ₀₀ . . . ϕ_(0N−1) . . . ϕ_(K−1N−1)) andwherein g1_(kn), g2_(kn) and g3_(kn) are multiplicative positiveamplification functions.
 36. The method of claim 34, wherein g1_(kn) isa positive non-linear amplification function with an amplification thatgradually decreases with increasing argument values.
 37. The method ofclaim 34, wherein g2_(kn) have a lower amplification for an argument of90 degrees than for an argument 0 degrees and that it is symmetricaround 90 degrees.
 38. The method of claim 34, whereing2_(kn)(ϕ_(kn)(i,j)) is limited to a band of values defined between anupper and lower limit.
 39. The method of claim 34, wherein g3_(kn) is afunction fitted to a training set and with (A₀₀ . . . A_(0N−1) . . .A_(K−1N−1), ϕ₀₀ . . . ϕ_(0N−1) . . . ϕ_(K−1N−1)) as argument.
 40. Themethod of claim 34, wherein g3_(kn)(A₀₀ . . . A_(0N−1) . . . A_(K−1N−1),ϕ₀₀ . . . ϕ_(0N−1) . . . ϕ_(K−1N−1))=g3_(kn)(A_(k0) . . . A_(kN−1)) andwherein amplification g3_(kn) is larger for orientations n close ton_(max), wherein n_(max) is calculated as the orientation with maximumamplitude.
 41. The method of claim 25, wherein a parameter a controlsthe variation of g3_(kn) for different n.
 42. The method of claim 41,wherein α is a function of grayscale, segmentation value, the spread ofn_(max) in a M×M neighbourhood or any combination.
 43. The method ofclaim 25, wherein said high-pass filtered detail image h₀ is non-linearprocessed by: (a) calculating for said pixel, at least one statisticalmeasure for two or more translation difference image pixel values withina neighbourhood of said pixel; and (b) modifying the value of said pixelof said detail image h′₀ as a function of said statistical measure,amplitude and phase of said detail image b_(kn) and value of said pixelof said detail image h₀, as follows:h′ ₀(i,j)=a ₀(i,j)h ₀(i,j)wherein,a ₀(i,j)=ƒ₀(√{square root over (var₀)},A _(kn),ϕ_(kn)).
 44. The methodof claim 25, wherein said high-pass filtered detail image h₀ isnon-linear processed byg ₀(i,j)=h ₀ *g ₀(A ₀₀ . . . A _(0N−1) . . . A _(K−1N−1),ϕ₀₀ . . .ϕ_(0N−1) . . . ϕ_(K−1N−1)).
 45. The method of claim 25, whereinprocessing at least one pixel of said detail images is repeated at leastone time for already processed detail image pixels.